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Supplementary Information 


Reference baseline: Central to the planetary boundaries framework is a comparison of the 
state of the individual boundary processes to Holocene conditions. In referring to the 
Holocene, we use its geological definition, as the latest warm inter-glacial Epoch of the 
Quaternary Period (covering the last 3 million years) of the Earth system. Technically, the 
Holocene covers the last 11,700 years of Earth's history and extends to the present. For the 
purposes of this study, however, we refer to Holocene conditions as those the Earth 
experienced from 11,700 years ago until 0 BCE or 1700, depending on the analysis. This pre- 
industrial period in Earth history constituted an extraordinarily stable state of the Earth 
system, with a climate system oscillating only 0.5 °C from the pre-industrial global mean 
surface temperature (9). During the whole of the Pleistocene epoch (from 2.6 million years 
ago), the Earth system has oscillated between long glacial periods followed by shorter inter- 
glacial intervals, with a number of inter-glacial periods identified in just the last one million 
years, of which the last is the Holocene) (/9, 98). Earth system modelling and ice-core data 
show that the Earth system never reached a change of 2°C of global mean surface 
temperature during all these warm inter-glacials (99), providing evidence of the Earth 
system's capacity to remain within inter-glacial conditions, to a certain degree even outside of 
the Holocene range, if Earth system interactions and feedbacks are kept intact. "Holocene- 
like" is then used in the text to refer to an Earth system that may be outside of the Holocene 
maximum range of variability for Earth system regulating processes and systems (planetary 
boundaries) but still remains in an inter-glacial state, i.e., has not crossed thresholds that 
irreversibly move the Earth system away from Pleistocene inter-glacial conditions (e.g., > 
2°C increase of global mean surface temperature). The reference point for a safe operating 
space of the Planetary Boundaries, is thus to maintain the Earth system in a manageable inter- 
glacial Holocene-like state. 


Boundary updates: 
Biosphere Integrity: 


Genetic diversity within species facilitates adaptation and evolutionary change under 
changing environments (/00,/0/) and promotes the resilience, productivity, and functionality 
of both marine and terrestrial ecosystems (/02, 103). However, implementing a control 
variable for genetic diversity in the Planetary Boundaries framework has been hampered both 
by lack of data, adequate predictive methods and of a consensus on what are the most 
relevant genetic metrics. The recent accumulation of genetic data across thousands of species 
from the poles to the tropics, the revolution of whole-genomic sequencing and the inception 
of new approaches, i.e., macrogenetics (104-108) pave the road for capturing the global 
genetic variability and its changes over time. Currently, knowledge of the global geography 
and the temporal trends of intra-specific diversity is still limited for many taxa across the tree 
of life (104, 106) and methods to predict future changes are still in development. Ideally, 
however, within a decade we should be in the position to estimate across hundreds of species 
losses in haplotypes, heterozygosity, functional genomic diversity or mutation load and 
predicting future changes using Mutation-Area Relationships MARs (23) for capturing 
globally extinction dynamics below species level (/09) and therefore surpass the inability of 
species-level extinction estimates to capture on-going rates of biodiversity loss. 


Biosphere function: Functional integrity of the biosphere is provisionally (i.e., until a more 
comprehensive metric of the functional influence of the biosphere on the state of the planet is 
available) measured by the human appropriation of net primary production, HANPP (cf. 32) 
in comparison to the Holocene value and variation of NPP. NPP is used as a proxy for the 
energy available to the biosphere to support its growth, maintenance, differentiation, and 
networking. HANPP includes both NPP altered (largely: reduced or eliminated) by human 
activity or harvested (minus backflows). 


The value and the variation of NPP during the pre-industrial Holocene was computed as 
follows. A 130,000-year climate simulation of the last glacial cycle, including the Holocene, 
was conducted with Climber-2, an Earth system model of intermediate complexity (//0). The 
time step of the simulation was 100 years and atmospheric COz concentration prescribed. 
This climate data set was applied to drive the biogeochemical global dynamic vegetation 
model LPJmL4 (95,///) to compute the evolution of terrestrial biome distributions, carbon 
and water cycle. Values of NPP were analyzed for the period 10,000 to 0 BCE, covering the 
bulk of the stable pre-industrial Holocene, by statistically combining the 100-year climate 
inputs with annual variation taken from historical CRU TS 3.23 climatology (//2,113), as 
realistic vegetation simulation requires natural variability of climate. For evaluation, 
however, the resulting annual values of Holocene NPP were subsequently once more 
averaged for each century to produce a mean NPP on the time scale of a century, removing 
the annual variation of climate applied. A histogram of the resulting century-scale Holocene 
NPPs was derived as well as the mean, median and standard deviation (the histogram 
approximates a bell-shape form, Fig. $1). 


Performance of the Climber-2 climate simulation for the Holocene was evaluated by 
comparing the histogram of temperature variations with the reconstructions (9, //4). While 
Kaufman et al. (74) combine both marine and terrestrial temperature proxies, albeit with 
some gaps in geographic coverage, Osman et al. (9) combine only marine proxies with the 
dynamic properties of a climate model to derive estimates of climate. While the Kaufman 
study (//4) produces slightly declining global mean temperatures after the middle of the 
Holocene, the study by Osman et al. (9) shows quite constant temperatures with no decline in 
that period. These authors attribute the decline seen in Kaufman et al. (//4) and previous 
estimates to the method used for deriving global temperatures from the proxies and not their 
climate model. The Climber-2 simulations agree in the evolution of temperature very well 
with the findings of Osman et al. (9), showing a climatically very stable Holocene after the 
initial transition phase. The temperature histograms show that the Climber-2 simulations and 
both reconstructions have a very similar variation of temperatures occurring over the 
Holocene. Note, however, that the reference period for these data sets is different; therefore, 
the absolute values of the histogram means cannot be directly compared. The figure shows 
these histograms (left) and the histogram of NPPs from the Climber-2 simulation as applied 
to the LPJmL biosphere model (right), also showing the relationship between histograms of 
annual NPPs (but with modern variability), and the century-mean NPPs used for subsequent 
analysis. Holocene mean NPP was estimated from the modelling to have been 55.9 GtC yr", 
the median 55.8 GtC yr", and the formal 1o-standard deviation 0.54 GtC yr’. 


Historical simulations of NPP: The LPJmL4 biosphere model is a well-established, 
documented, and evaluated dynamic global model of terrestrial biogeochemistry and 
vegetation structure that is widely used. It allows simulation of vegetation and soil carbon 
and water fluxes for both natural and fractional agricultural areas at a spatial resolution of 0.5 
degrees, amounting to nearly 60,000 pixels to cover the relevant global land surface. See 


Schaphoff et al. studies (95, 1/1) for details of model structure, process representations and 
benchmarking. 


Daily computations are processed to generate annual output for the simulations used here. A 
20,000-year spin-up for the carbon pools produced by natural vegetation is followed by the 
130,000-year simulation of the last glacial cycle; a 3000-year spin-up is followed for the 
historical simulations by a spin-up introducing land use patterns from 1500. The simulation 
then follows dynamic expansion of land use from 1700 to the present using concurrent values 
of atmospheric CO2 concentration. However, climate for the years 1700-1900 is a statistical 
sample of the period 1901-1930. After 1901, the simulation uses observed CRU climatology, 
actual land use and atmospheric CO2 concentrations until 2009. Evaluations are therefore 
restricted to the time from 1900, except for an estimate of NPP for 1700. Simulations were 
performed for potential natural vegetation (without any land use), for actual land use, and for 
keeping atmospheric CO2 concentration constant at the pre-industrial value in order to 
separate the direct effects of CO2 fertilization. Resulting annual NPP was averaged over a 30- 
year period and the resulting smooth evolution fitted for interpolation and producing an 
estimate for the year 2020. 


HANPP: Quantitative evaluations of human appropriation of NPP (HANPP) are scarce in the 
literature and vary due to differences in definition, detail, and data used. A review of 
approaches and findings is given in Haber] et al. (32). Fig. S2 illustrates the authoritative 
definition of HANPP used here. The most detailed and systematic derivation is found in 
Haberl et al. (1/5) for the year 2000, and a subsequent evaluation for the 20" century in 
Krausmann et al. (33). Although these were derived from baseline NPP simulations for 
potential natural vegetation also produced by the LPJmL model (from which human 
appropriation was deduced using a variety of data bases of human land use in a number of 
relevant categories, corrected for backflows and fires), the model version used then was not 
the same as the more advanced version of LPJmL used here. We therefore applied the values 
of HANPP as a percentage of potential natural NPP (33) to the potential natural NPP values 
of the LPJmL simulations used here to derive absolute values of HANPP for selected years. 
As these displayed a very smooth evolution, they were fitted to derive intermediate values 
and an estimate for the year 2020. 


Evaluation: Fig. S3 and Table S1 present an overview of key NPP results from the biosphere 
model simulations and corresponding HANPP values. The position of the (provisional) 
planetary boundary component for functional Biosphere Integrity is marked at 10% and the 
upper end of the zone of increasing risk at 20% of the pre-industrial Holocene mean value of 
terrestrial NPP. 


Climate Change: 


Climate Change control variables of atmospheric CO2 concentration and total effective 
anthropogenic radiative forcing and their boundary levels are retained (/,2). The planetary 
boundary for atmospheric CO2 concentration is set at 350 ppm and for radiative forcing at | 
W m”. The upper end of the zone of increasing risk is set at 450 ppm for atmospheric CO2 
concentration and at 1.5 W m” for radiative forcing (2). The control variable of radiative 
forcing describes the aggregate most important anthropogenic impact on the Earth’s energy 
balance via release of greenhouse gases, aerosols, and albedo forcing (17), while the control 
variable of atmospheric CO2 concentration highlights the currently dominant role and policy 
relevance of this greenhouse gas. The boundary value for total effective anthropogenic 


radiative forcing of 1 W m” is chosen such that it corresponds to the radiative forcing 
contribution of atmospheric COz alone at its boundary value of 350 ppm (that was reached 
between 1987—1988), when the positive and negative contributions of other greenhouse 
gases, aerosols and land-use changes had approximately cancelled out (/7). 


Ozone Depletion 


Transgression and current ozone values were estimated from daily measurements from 
NOAA (96) for both hemispheres. These data were spatially aggregated into global monthly 
medians. In order to compensate for the high annual fluctuations, the current value was 
calculated as a five-year median for 2016 - 2021. 


Freshwater Change 


Freshwater supports terrestrial and aquatic ecosystems, mediates flows of material and 
energy, and regulates the Earth’s climate. Thus, it is tightly coupled with the core boundaries 
of Climate Change, Biosphere Integrity and Novel Entities, as well as the Land-System 
Change and Biogeochemical Flows boundaries (14,22). While the PBs for Novel Entities and 
Biogeochemical Flows implicitly incorporate quality aspects of freshwater, the Freshwater 
Change PB explicitly addresses the quantity of freshwater. Terrestrial and freshwater 
ecosystems have generally adapted to specific quantities of freshwater, which vary naturally 
both within and between years. The level of wetness of landscapes also regulates climate 
from micro- to regional and global scales, such that changes in wetness of one location can 
impact rainfall and consequently river flows and soil moisture both locally and in distant 
locations (J 16-118). 


The earlier PB for Freshwater use explicitly accounted only for human consumptive use of 
blue water (i.e., from rivers, lakes, reservoirs, and renewable groundwater stores) as a proxy 
for water cycle change over land, with minimum levels of streamflow requirements in aquatic 
ecosystems as a sub-global boundary condition (2,//9). This control variable misses — or 
aggregates too broadly — various important anthropogenic disruptions to the hydro-ecological 
and hydro-climatic Earth system functions of freshwater, such as land-use, land management 
and climate change impacts on soil moisture, evaporation and precipitation (/ 16,120,121), 
and dam operation impacts on temporal variability of streamflow (/22). The new approach 
chosen here better acknowledges the many Earth system impacts that freshwater changes 
have beyond aquatic ecosystem integrity on, for example, terrestrial productivity, biome 
resilience (/23), carbon and nutrient transport, uptake, and storage (124,125), as well as 
regional atmospheric circulation (/26) by accounting for changes in both green water (soil 
moisture) and blue water (streamflow). Both local and global limits for freshwater change are 
based on variability within the anthropogenically relatively undisturbed pre-industrial period 
(used as proxy for Holocene conditions), such that persistent and widespread exceedance of 
these bounds represents a shift in the freshwater cycle. 


Boundary value and current states of the Freshwater Change boundary were derived by 
analyzing spatially and temporally explicit streamflow (blue water) and root-zone soil 
moisture (green water) data from the inter-model comparison project ISIMIP 2b (127). For 
root-zone soil moisture, data from four global hydrological models (GHMs) were used: 
CLMS0, LPJmL, MPI-HM, and PCR-GLOBWB. Streamflow data were derived from six 
GHMs: H08, LPJmL, MATSIRO, MPI-HM, PCR-GLOBWB and WaterGAP2. Each GHM 
was forced with modelled climate from four different general circulation models: GFDL- 


ESM2M, HadGEM?-ES, IPSL-CM5A-LR and MIROCS; with the exception of MPI-HM 
lacking HadGEM2-ES. This resulted in 15 ensemble members for root-zone soil moisture 
and 23 for streamflow. 


The control variables were calculated as follows, for green and blue water separately 

1. The local baseline range was set by determining dry (Sth percentile) and wet (95th 
percentile) bounds during the pre-industrial (1661—1860) period, for each month of 
the year, for each 30x30 arcmin grid cell and for each ensemble member. Values 
outside of these bounds are referred to as “local deviations”. 

2. For each month and each ensemble member, the percentage of global ice-free land 
area with local blue/green water deviations, both for dry and wet deviations, was 
calculated by aggregating the areas of grid cells with local deviations. Areas with dry 
and wet deviations were summed to derive total area with deviations. 

3. For each year and ensemble member, the mean of monthly areas with deviations was 
used to calculate the percentage of annual global ice-free land area with local 
blue/green water deviations. 

4. For each year, the median of ensemble members’ values was taken. 


These analyses result in a time series of annual global land area with streamflow/root-zone 
soil moisture deviations, starting from the pre-industrial period (1661—1860) and continuing 
to the industrial period (1861—2005). The end date is limited by the ISIMIP 2b simulation 
protocol, due to lack of reliable estimates on some of the major anthropogenic forcings, such 
as irrigation extent (/28), at the appropriate scales. Global trends in many of the key drivers, 
such as irrigation area (84), water use (/29), dam construction (/30) and forest loss (131) 
have increased since 2005. Consequently, the “current’’ estimate (mean of 1996—2005) 
presented here is likely conservative. 


The annual global land area with local deviations during the pre-industrial period is used to 
define the global boundary condition. The global land area with local blue and green water 
deviations was relatively stable during this period, varying approximately within 1.5 
percentage points around the median of 9.4% (blue) and 9.8% (green). The global boundary 
is provisionally set at the upper limit (95" percentile) of this pre-industrial variability, which 
was 10.2% for blue water and 11.1% for green water. This precautionary placement was 
chosen because scientific evidence of the levels of blue and green water deviations that may 
trigger irreversible changes in Earth system functions is currently lacking. To acknowledge 
uncertainties related to both data and the exact boundary position, the upper limit of the zone 
of increasing risk is provisionally set at 50% of land area with freshwater deviations, for both 
blue and green water, following the reasoning by Barnosky et al. (/32). At the landscape 
scale, undisturbed patches are prone to tip into a new state when 50-90% of surrounding 
patches are disturbed. Given the interacting and multiplying effects at larger scales, shifts in 
50% (or less) of the Earth’s ecosystems can be expected to trigger global forcings and rapid 
changes in remaining systems (/32). 


Currently (mean of years 1996-2005), about 18.2% (blue) and 15.8% (green) of the global 
land area experience local freshwater deviations, which we interpret as a substantial 
exceedance of the Freshwater Change boundary, given the very narrow variability in the pre- 
industrial period. The blue water boundary was transgressed in 1905, and the green water 
boundary in 1929 (46). Even the ensemble interquartile ranges surpass the boundaries in the 
latter half of the 20" century, which increases confidence on our conclusion that the boundary 
is transgressed but also underlines the uncertainty on the exact timing and magnitude of 


transgression (46). Transgressing the Freshwater Change boundary is consistent with ‘the 
great acceleration’ (2) of anthropogenic drivers impacting freshwater. Our estimates of an 
increase in freshwater deviations are in good agreement with results from many other studies 
of changes in the global water cycle. These include increases in precipitation (/33) and river 
flow (134) extremes, increases in drought occurrence and severity (/35), intensification of the 
water cycle (136), widespread and severe violations of environmental flows (136-139) of 
river flow regimes (/22). The geographic patterns identified in our assessment also agree well 
with previously reported research, where tropics and subtropics show the greatest drying and 
the northern latitudes the greatest wetting (46, 140). 


The conclusion that human activities have brought the Freshwater Change planetary 
boundary into a zone of increasing risk of changing Earth system state is consistent with 
many reported freshwater-mediated impacts. For blue water, perhaps the most dramatic 
example of a link to Earth system processes is the destruction of the Aral Sea, which occurred 
due to substantial streamflow reductions from overuse of water for irrigation, and resulted in 
lake depletion with consequent ecological degradation and regional climate change (/4/). 
Green water changes have been associated with productivity loss in both natural and 
cultivated lands, as exemplified by drying-induced forest dieback in all continents (/42) and 
increases in flood and drought-induced food production shocks, particularly in South and 
East Asia, Australia, and North Africa (143). Examples of ecological and climatic impacts of 
water surpluses are habitat loss in the Central Amazon floodplains due to anthropogenic flood 
pulse disturbance (/44), and increased greenhouse gas emissions from reservoirs (145), both 
created by dam construction and operation. Reversing the Freshwater change PB 
transgressions will likely require reduced human disturbance through land-use and land 
management, water withdrawals, large water transfers, reservoirs and dam operations, and 
greenhouse gas and aerosol emissions. Note, however, if the previously defined boundary 
(/,2) were applied, human impacts on this boundary would still be considered as being within 
the safe operating space (49). 


Aerosol Loading 


The atmospheric burden of small particles (aerosols) have both physical and biogeochemical 
effects in the Earth system. They influence Earth’s energy balance by absorbing and 
scattering of radiation, and changing cloud properties (/7). These physical effects are 
represented in the climate change planetary boundary, but there are also other complex 
interactions between aerosols, climate and the biosphere. One example is the dimming effect 
of aerosols, i.e., a higher aerosol loading leads to less surface irradiance, which can cause for 
instance a relative reduction in surface warming and evaporation. Aerosols further affect the 
hydrological cycle by influencing clouds and precipitation (/46), although these are yet to be 
fully understood. Moreover, aerosols influence primary production, human and organism 
health, and ecosystem function, e.g., through the deposition or inhalation of trace substances, 
and altering the direct irradiance that affects photosynthesis. 


Aerosol information is here based on the data assessment by (57). We calculate the area- 


weighted hemispheric difference (AT) in the aerosol optical depth (t) for present-day. T is a 
measure of the aerosol extinction, i.e., scattering and absorption of irradiance. For being close 


to observations, we use all twelve t datasets from satellites, re-analysis, and climatologies in 
the time period 1998-2020 that are listed in Table 1 of (57). We first calculate for each data 

set At = t(@) — tT(—@) at each latitude (@) between 60 degrees North and 60 degrees South. 
We than compute the average over all latitudes weighted by the latitude-dependent area. The 


reported values for AT are the mean across all datasets and the data-to-data standard deviation 


as uncertainty measure. We assume that aerosol impacts scale with At. Given that monsoon 
dynamics and the associated regional rainfalls respond to changes in anthropogenic aerosols 
(147-150), we assume the effects of anthropogenic aerosol release mimic reseale from natural 
sources. 


Ocean Acidification 


The current (i.e., year 2022) global averaged saturation state of aragonite (Q=2.8) is 
estimated from the climatological average value of 3.03 in year 2000 and the corresponding 
global decrease of 0.1 per decade (7/). This corresponds to an estimated decrease of 0.64, or 
81%, of the pre-industrial value of 3.44 (2). 


Land System Change 

This update retains the same general approach used earlier (2). However, improved 
methodology is available for estimating the remaining area of forest. These improvements 
include major advancements in satellite image technology and improved classification 
algorithms, which formerly relied on statistical methods but have now shifted towards using 
machine learning approaches. 


Global: The global boundary is defined as the area of forested land (on the ice-free land 
surface) that is maintained, expressed as a percentage of the potential area of forested land in 
the Holocene (that is, the area of forest remaining assuming no human land-cover change). 
The boundary of 75% of potential forest cover remaining is retained. The boundary has been 
constructed as a weighted aggregate of the three individual biome boundaries as described 
below. 


Biome: Boundaries for each forest biome are defined as the forest area maintained in each of 
the three major forest biomes — tropical, temperate, boreal — expressed as a percentage of the 
potential forest area in each of these three biomes. The estimated boundary for each of the 
biomes is based on (i) the relative potential of land-cover change within each biome to 
influence the climate system at the global level, and for some biomes, (ii) the potential for a 
threshold of land-cover change beyond which self-reinforcing feedbacks within the biome 
lead to land-cover change across a much larger area. This approach is identical to the one 
used in 2015. 


The boundary levels and current status of the eight major forest biomes are given in Table 1. 
Seven of the eight biome-level boundaries have now been transgressed, one more than in the 
2015 assessment. In addition, there have been some notable shifts in the amounts and patterns 
of deforestation: 


Tropical forest. The boundary is retained at 85% of potential forest cover to be maintained, 
based on the potential for a threshold of land-cover change in the Amazon Basin tropical 
forest beyond which much of the forest is converted to a savanna or grassland. The 
mechanism for such a transition, described in more detail (2) but is based on the reduction of 


evapotranspiration as the tropical forest is converted into grazing or cropping land, and the 
subsequent decrease in rainfall across the basin. 


Temperate Forest. There has been no change in the overall pattern of deforestation of 
temperate forests compared to the 2015 assessment (2). Both European and Asian forests 
remain in a transgressed state, while North American temperate forests are still within the 
safe operating space. However, the areas of both the North American and European 
temperate forests appear to be considerably reduced. The Asian temperate forest, on the other 
hand, has increased in area from since 2015. While changes in the methodology used to 
estimate forest cover may explain some of these differences, the some of the increase in 
forest cover in Asia may be due to the very large reforestation initiative in China. 


Boreal forest. Both the North American and Eurasian boreal forests have decreased very 
slightly in area since 2015. Given the changes in methodology, these estimated increases in 
boreal forest cover are not likely to be significant. The current areas of both boreal forests 
indicate substantial transgression of this boundary. 


Earth system effects of differing scenarios of transgression of Land System Change and 
Climate Boundaries 


We apply ESM CM2Mc-LPJmL (a configuration of the Potsdam Earth Model POEM (85), 
which includes marine (BLING) and terrestrial (LPJmL5) biosphere models, to explore short- 
and long-term biosphere responses to climate change. Especially the model’s advanced 
terrestrial biosphere is able to account for important biological processes such as competition 
between vegetation functional types, tissue mortality, fire dynamics and water stress that 
influence the land carbon balance and biophysical feedbacks on climate. 


Development of terrestrial biologically mediated C-sinks: 


CM2Mc-LPJmL v1.0 (85) is a coupled Earth system model, combining the coarse-grained 
but relatively fast atmosphere and ocean model CM2Mc (/5/) with the state-of-the-art 
dynamic global vegetation model (DGVM) LPJmLS5 (95, //1,, 152). CM2Mc is a coarse 
(3°x3.75° latitude-longitude) configuration of the Climate Model CM2 (/53), developed at 
the Geophysical Fluid Dynamics Laboratory (GFDL). The original CM2Mc model includes 
the Modular Ocean Model 5 (MOMS) and the global atmosphere and land models AM2-LM2 
or AM2-LM (/54) with static vegetation. In CM2Mc-LPJmL, the land component LM/LM2 
is replaced by the dynamic global vegetation model LPJmL5. AM2 and MOMS remains 
dynamically coupled to the model framework. The different model components are connected 
by the Flexible Modeling System (FMS), developed by GFDL. 


LPJmLS5 (Lund-Potsdam-Jena managed land) is a state-of-the-art DGVM, which is well 
established and thoroughly validated (95,//1,152)). It models water and carbon fluxes for 
natural and managed land, carbon stocks and the global surface energy balance, but the 
nitrogen cycle remained deactivated in keeping with previous model setups ($5). LPJmL5 
simulates the impact of bioclimatic limits and effects of heat, productivity, and fire on plant 
mortality to determine the establishment, growth, competition, and mortality for different 
plant functional types (PFTs) in natural vegetation and crop functional types (CFTs) on 
managed land. The model is forced by climate and soil data. New developments since the 
original publication (/55) led to the incorporation of a water balance (/56), agriculture (157), 


wildfire in natural vegetation (/58,/59), and the impact of multiple climate drivers on 
phenology (/60,/6/). 


Using LPJmLS, enables the model to make, e.g., use of its advanced land use scheme, the 
process-based fire model SPITFIRE, a representation of permafrost and a state-of-the-art 
water cycling. These processes within LPJmLS have large effects on the location, timing and 
magnitude of the atmospheric water fluxes and precipitation through plant 
evapotranspiration, surface temperature and the simulated canopy humidity. The fully 
coupled energy and water cycle allows to investigate the impact of biophysical atmosphere- 
biosphere feedbacks on global climate trajectories and to quantify the impacts of 
deforestation or afforestation scenarios. 


Compared to the published version of CM2Mc-LPJmLv1.0 (85), the version used in this 
paper incorporates a few small changes, including setting the net water influx from rivers into 
the ocean to zero to avoid unrealistic high fluctuations of the sea level. The sea level can still 
change due to sea ice melt and thermal expansion. We also improved the initialization of the 
model after starting from a previous state by including more variables in the saved state file. 
Finally, we slightly adapted the canopy conductance used in the photosynthesis routine to 
increase plant productivity and better match historic trajectories of the land carbon sink. 


Model protocol: 


The model experiments of this paper are consistent with Driike et al. (85): After a 5000-year 
stand-alone LPJmL spin-up, a second fully coupled spin-up under pre-industrial conditions 
without land use was performed for 1500 model years. In this way, we ensure that the model 
starts from a consistent equilibrium between the long-term soil carbon pool, vegetation, 
ocean, and climate. 


The subsequent transient historic phase of the model is performed from 1700-2018, using 
historic land use data from 1700 (162) and historic concentrations of greenhouse gases, solar 
radiation, ozone concentrations and aerosols from 1860, which were kept at pre-industrial 
conditions beforehand. From 2004 on, only greenhouse gas forcing remains, while aerosols, 
solar radiation and ozone are set to their corresponding 2003 values. Scenario values for 
atmospheric CO2 concentration were set to 350, 450 and 550 ppm, while scenarios for land 
system change (pattern of remaining forest) were generated by an algorithm that replaces 
natural forest by managed grassland and crops adjacent to land-use areas until the extent of 
natural land was reduced to the designated fractions in the different climate zones (tropical, 
temperate, boreal). 


In the following scenario names, climate change is denoted CC and land system change as 
LSC. The descriptors “PB”, “upper” and “beyond” refer to three levels of forcing of these 
planetary boundary dimensions. The “planetary boundaries (PB)” setting delineates the safe 
operating space, i.e. the planetary boundary itself; the “upper” value is at the upper end of the 
zone of increasing risk, identified based on conditions that while beyond the boundary are 
still potentially in an interglacial state; and the highest value analyzed (“beyond”) gives a 
forcing that likely endangers the interglacial condition of Earth. 


e PBCC and PB LSC: From 1989 on constant atmospheric CO2 concentrations at 350 
ppm and constant land use pattern from 1989 for 800 years. 


e Upper CC and PB LSC: From 2019 for 10 years linear trajectory to 85/50/85 % 
(tropical/temperate/boreal) remaining natural vegetation and to 450 ppm (1 % yearly 
increase from 2019 levels). From 2029, constant land use pattern and constant 
atmospheric CO2 concentration at 450 ppm for 800 years. 

e Upper CC and upper LSC: From 2019 for 10 years, linear trajectory to 60/30/60 % 
(tropical/temperate/boreal) remaining natural vegetation and to 450 ppm (1 % yearly 
increase from 2019 levels). From 2029, constant land use pattern and constant 
atmospheric CO2 concentration at 450 ppm for 800 years. 

e Upper CC and beyond LSC: From 2019 for 10 years, linear trajectory to 40/20/40 % 
(tropical/temperate/boreal) remaining natural vegetation and to 450 ppm (1 % yearly 
increase from 2019 levels). From 2029 constant land use pattern and constant 
atmospheric CO2 concentration at 450 ppm for 800 years. 

e Beyond CC and upper LSC: From 2019 for 34 years, linear trajectory to 60/30/60 % 
(tropical/temperate/boreal) remaining natural vegetation and to 550 ppm (1 % yearly 
increase from 2019 levels). From 2053, constant land use pattern and constant 
atmospheric CO2 concentration at 550 ppm for 800 years. 

e Beyond CC and beyond LSC: From 2019 for 34 years, linear trajectory to 40/20/40 % 
(tropical/temperate/boreal) remaining natural vegetation and to 550 ppm (1 % yearly 
increase from 2019 levels). From 2053, constant land use pattern and constant 
atmospheric CO: concentration at 550 ppm for 800 years. 


Model Evaluation: 


The performance (Fig. S6) of the newly coupled CM2Mc-LPJmL model is similar to 
CM2Mc-LaD (140) and comparable to results from the Coupled Model Intercomparison 
Project (CMIP5, CMIP6) (/63). Furthermore, the model shows a stable performance over 
750 years and reasonable reactions to climate and land-use change (for details see (S5)). 


Fig. S4 shows the modeled spatial patterns of al0-year average of recent historical surface 
temperature and precipition, and a comparison of their zonal means with ERAS reanalysis. 
Fig. SS shows the modeled temporal development of the yearly and decadal global mean 
temperature anomaly (relative to the reference period 1951—1980) compared to GISTEMP 
data. Fig. S6 compares the simulated distribution of above-ground biomass against the 
GlobBiomass data set (/64). 


Model results (Fig. S7, Fig. S8, Fig. S9, Table S2) 


The “safe” climate and land system change boundaries were passed around 1990. 
Atmospheric CO2 concentration passed 350 ppm and the remaining fraction of natural land 
area fell below 85/50/85 % of tropical/temperate/boreal forest cover (lower LSC). Today, we 
are heading to the upper end of the zone of increasing risk for these two boundaries. The 
“next stop” before entering dangerous territory is therefore the upper end of the zone of 
increasing risk at 450 ppm and the upper end of the zone of increasing risk for land system 
change with 60/30/60 % forest cover retained (upper LSC). 


To complete the modeling experiments, we also ran experiments maintaining the climate 
dimension at the upper end of the zone of increasing risk at 450 ppm, while the LSC 
boundary was maintained (LSC at the planetary boundary) or transgressed to beyond the zone 
(40/20/40 %, transgressed LSC), and experiments where CC was beyond the zone of 


increasing risk at 550 ppm and LSC was either at the upper limit of the zone of increasing 
risk or beyond. 


The modeling results provide support for the placement of the Climate Change planetary 
boundary at 350 ppm, as the temperature increase from pre-industrial would be within ca. 1.3 
°C above preindustrial levels, below the Paris target of 1.5 °C and only a small cumulative 
net amount of carbon (ca. 60 GtC) emitted from the biosphere over the modeling period of 
800 years. 


At the same time, the long-term change in terrestrial carbon stocks indicates that respecting 
the planetary boundary for LSC retains a stabilizing influence, i.e., a C sink, on climate when 
atmospheric COz is 450 ppm, in contrast, an LSC value advanced to the upper edge of the 
zone of increasing risk produces the opposite, an additional leakage of carbon to the 
atmosphere, i.e., an acceleration of warming and hence implies further transgression of 
climate into the dangerous zone. This is even more the case if the upper end of the zone if 
increasing risk for the land system change boundary is transgressed. 


If atmospheric CO2 concentration were to rise to 550 ppm, i.e., attain a value beyond the zone 
of increased risk for the Climate Change boundary, not only would an 0.7-0.8 °C of 
additional long-term warming occur over land compared to the long-term temperature at 450 
ppm, but also around 150 GtC of carbon loss would occur from land. The results of the 
different experiments are summarized in Table S2 and Fig. 2 of the main text. 


When transgressing the LSC boundary by expanding managed land in CM2Mc-LPJmL, 
natural vegetation is partially replaced by pasture and crops. The decrease of biomass and 
forest cover affect biophysical feedbacks mainly in three distinct ways: 1) Crops transpire 
less water than large forests, thus decreasing the water flux to the atmosphere and in turn 
humidity, precipitation, and cooling by latent heat exchange (/65). 2) Crops and grasses have 
a higher albedo than a closed forest, therefore impacting the energy cycles towards decreased 
surface temperature (/66). 3) Temperatures increase due to a decrease in roughness length 
(167). The net effect of these biophysical feedback mechanisms is a warmer and drier climate 
in most regions. Specifically in the tropics, latent and sensible heat fluxes far exceeded the 
potential cooling of an increased surface albedo. 


Development of biologically mediated C-sinks in the ocean: Model simulations of carbon 
fluxes, primary production (PP) and the saturation state of aragonite (Q) were analyzed in the 


initial (1.e., from 1988) and final (after ~800 years) 30-year period. In addition, biogenic 
particulate fluxes were compared with empirically derived estimates calculated from 
temperature and PP (Table S3). The biological module was based on the BLING model and 
monthly averaged fields of temperature, biological phosphorous uptake and the sinking flux 
of particulate phosphate (P) were applied for calculating PP and export fluxes of particulate 
organic carbon by assuming a molar C:P ratio of 106. New production (NP) and the 
mesopelagic C-flux (Fsoom) were calculated from the particulate flux at 100 m and 500 m 
depth, respectively. 


The empirically derived estimate of NP was calculated from sea surface temperature (SST), 
PP and the export ratio (ef) (1.e., ef = (a+b SST) PP*, where a,b and c are constants (90), and 
compared with simulated NP. Mesopelagic particulate C-flux (F) was estimated from the 
empirical relationship suggested by (9/) (i.e., F(z) = F(zo) (z/z0)°; F(zo) is the flux at 100 m 
depth, i.e., NP) where the power dependence of the remineralisation length scale (b) 


decreases linearly with temperature in the upper 500 m (To-soom). Monthly averaged values of 
the C-flux at 500 m depth (Fs00m-E) were calculated from NP and To-soom. 


Sensitivity of C-inventories and biogenic fluxes: SST increased by 1.4 °C between the 550 
and 350-simulation and temperature in the upper 500 m increased correspondingly by 1.1 °C 
(Table S3). The total inventory of DIC increased by 946 GtC of which 273 GtC accumulated 
in the upper 1000 m. The accumulation of DIC was mainly driven by increased atmospheric 
pCO? whereas only a minor fraction could be ascribed to biogenic fluxes. A characteristic 
response time (e-folding time) for accumulation of total DIC and DIC in the upper 1000 m 
was estimated (i.e., approximating three e-folding times with the elapsed time until the 
change in DIC had reached 95% of its final change) for the 550-simulation to be 261 and 155 
years, respectively. Both the 450 and 550 ppm simulations showed a continuous rise in DICtot 
during the entire 800-year period. 


New production decreased by ~4 % between the two experiments (550 and 350) whereas the 
decreased mesopelagic flux at 500 m depth of 0.08 and 0.19 GtC yr"! for the model and 
empirically derived fluxes, respectively, implied a significant reduction of 5-15% of the 
organic matter flux to the ocean interior. A small increase of NP between the 1988 and 350- 
experiment (in both the model and empirically derived NP) could be explained by the slight 
increase in PP (e.g., 51.3 GtC yr") in the end of the 350-experiment. Model- and empirical 
estimates of NP and Fsoom were generally in good accordance. 


Saturation state of aragonite: The initial Q in the model (1.e., calculated from the carbonate 
system and based on SST, surface salinity, DIC and alkalinity) was in general accordance 
with climatology, i.e., the simulated value of 3.2 in 1988 is within 0.05 of the estimated 
climatological value (Table S3). The final state in the 550 ppm simulation (Q=2.51) crosses 
the planetary boundary (80% of the pre-industrial value corresponding to Q=2.75). A 
characteristic response time was estimated (see above) and the 450 and 550-simulations were 
characterized by a decadal decrease of 0.1, in accordance with the current climatological 
change (7/). Thus, the major change in surface occurs within a few decades. 


Fig. $1. Histograms comparing the variation of temperature during the pre-industrial 
Holocene (10.000-0 BCE) for two reconstructions and the Climber-2 climate model 
simulation (left; note that the reference period for each period is different, so that the shape 
of the histograms can be compared but not the values of the anomalies plotted), and of the 
century-mean and annual (with modern climate variability) NPPs simulated with the LPJmL 
biosphere model (right). 
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Fig. S2. Schematic of the definition of HANPP used in this study (not to scale; following 
the definitions of (75)). 
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Fig. S3. Summary of key NPP and HANPP results from the LP JmL4 biosphere model 
simulations. 


Holocene Actual NPP w/o Potential NPP on Actual NPP _— Potential NPP 
NPP - HANPP CO, fertilization NPP-HANPP natural areas withland use w/o land use 
SS Holocene 


2020 
(extrapolated) 


I 
I 
I 
I 2000 
I 
I 
I 


\ 
\ 
1 
: \ 
\ 
20% 10% e bg 1700 
Percentage of (PB) 
Holocene NPP Mean 
Median 


Holocene century-mean NPPs 


38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70 72 74 


Terrestrial NPP (GtC/year) 


Fig. S4. (a) Global surface temperature from CM2Mc-LPJmL averaged over the period 
1994-2003, (b) global precipitation from CM2Mc-LPJmL averaged over the period 1994— 
2003, (c) zonal mean temperature from CM2Mc-LPJmL (red line) and ERAS data (blue line) 
averaged over the period 1994-2003, (d) zonal mean temperature from CM2Mc-LPJmL (red 
line) and ERAS data (blue line) averaged over the period 1994—2003. Plot is analogous to 
(85). 
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Fig. S5. Yearly and decadal global mean temperature anomaly (relative to the reference 
period 1951—1980) of CM2Mc-LPJmL compared to GISTEMP data from 1880-2018. Note 
that, from 2004 on, only greenhouse gas forcing remains, while aerosols, solar radiation and 
ozone are set to their corresponding 2003 values. Plot is analogous to (85). 
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Fig. S6. (a) Mean global above-ground biomass of GlobBiomass (198) evaluation data. (b) 
Mean global above-ground biomass of CM2Mc-LPJmL over the period 2006-2015. (c) 
Difference of the above-ground biomass between CM2Mc-LPJmL and GlobBiomass 
evaluation data. Blue/red colors denote an overestimation/underestimation of biomass by 
CM2Mc-LPJmL. (d) Latitudinal sum of above-ground biomass from CM2Mc-LPJmL (blue 
line, R2=0.65, NME=0.58) and GlobBiomass evaluation data (red line). Plot is analogous to 
(85). 
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Fig. S7. (a) Total land carbon, (b) vegetation carbon, (c) soil carbon and (d) litter carbon for 
the different experiments, including the historic period until 2019 in black. In the first 100 
years of the scenarios, the fertilization effect leads to a strong increase in total land carbon for 
the experiments at 450 ppm/lower LUC and 550 ppm/upper LUC, respectively. Soil carbon 
in the other experiments decreases due to increased heat stress and increased LUC. 
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Fig. S8. Global mean land temperature in a 30-year running mean for the different 
experiments, including the historic period until 2019 in black. Temperature increases mainly 
due to increasing atmospheric CO: concentrations and to a lesser extent due to biophysical 
feedbacks for the different land use change scenarios. Long-term temperature increase is due 
to lag effect of heat transfer in the ocean and to a lesser extent in the different soil layers. 
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Fig. S9. Global mean land precipitation in a 30-year running mean for the different 
experiments, including the historic period until 2019 in black. Precipitation increases mainly 
due to increasing temperatures, which increases evapotranspiration and humidity in the air. 
This effect is offset by a decrease of evapotranspiration in managed land. Therefore, a loss of 
natural land, leads to a decrease in evapotranspiration, humidity and consequently 
precipitation. 
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Table S1. Summary of key NPP and HANP?P results from the LPJmL4 biosphere model 
simulations. 


NPP (GtC/yr) PNV Actual on non- Actual with no HANPP PNV-HANPP Holocene 
(nolanduse) (withlanduse) agricultural co2 Mean - HANPP 
areas fertilisation 
Holocene 55.9 

1700 56.2 54.7 

1925 60.0 57.4 53.6 54.4 8.1 51.9 47.8 

1950 62.1 59.0 54.9 54.5 9.7 52.4 46.2 

1975 64.8 61.0 56.7 54.2 12.4 52.4 43.5 

2000 68.2 63.4 58.9 53.6 14.6 53.7 41.4 

2020 71.4 65.8 61.0 52.8 16.8 54.6 39.2 


Table S2: Land carbon exchange with the atmosphere (negative — land sink, positive — 
land source) and temperature increase over land and globally between the years 1988 and 
2100 (green) and 1988 and 2770 (orange). Note: Additional terrestrial carbon stock increase 
from pre-industrial to 1988 is 34.75 GtC, temperature warming from pre-industrial to 1988 is 
+0.58°C over land and 0.54°C globally. 


Scenario: CC / LSC CarboninGt CarboninGt Landtempin°C Land temp in °C Global temp in °C Global temp in °C 


350 ppm PB / PB +25.1 +68.3 
450 ppm upper / PB -127.4 -132.1 


450 ppm upper / upper +77 +80.8 
450 ppm upper / beyond +130.7 +210.5 
550 ppm beyond / upper -54.9 +6,1 

550 ppm beyond / beyond +63.2 +145.3 


Table S3. Model simulations and empirically derived results (initially (1988-2018) and 
at the end of the three experiments): Sea surface temperature (SST), average temperature in 
the upper 500 m (To-soom), primary production (PP), new production (NP), particulate organic 
C-flux at 500 m depth (Fsoom), total ocean inventory of DIC (DICiot) and in the upper 1000 m 
(DICio00m), and the saturation state of aragonite (Q). Empirically derived fluxes (NP-E, Fsoom- 
E) are shown for new production and particulate organic C-export at 500 m depth. 


Exp. SST To-500m PP NP NP-E  Fs00m ~~ Fso0m-E DICtor DIC i1000mQ 


ppm °C “C GtC yr! GtC yr! GtC yr! GtC yr! GtC yr'EgC Eg C 


ini(1988)18.8 11.3 S302 8.5 8.8 1.59 170 384 9.1 3.2 


350 19.1 11.7 51.3 8.6 9.0 1.62 173 387 9.1 3.2 
450 19.8 12.30 51.3 835 8.9 1,59 164 39.1 9.2 2.8 


550 20.5 12.8 51.3 83 8.7 1.54 [54 .393. 93 Zz 
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